
*Import the optimpal policy solutions: 
use "$root/Results/estimation_results/opt_pol_byquart.dta", clear

**Add the initial calculations to assess the changes in the outcomes of interest: 
frame create init
frame change init
import excel "$root/Results/estimation_results/counterfactual/init_prob_byinc.xlsx", firstrow clear

gen prob_EV=prob if drivetype==3


gcollapse (sum) co2_ann emb_co2_ann petrol_tax veh_tax tax_total tax_pv prob_EV, by(wealth_group)

expand 2 if _n==1, gen(newvar)

replace wealth_group=0 if newvar==1
drop newvar


foreach var of varlist co2_ann emb_co2_ann petrol_tax veh_tax tax_total tax_pv {
	
	replace `var'=. if wealth_group==0
	egen tot_`var'=sum(`var') 
	replace `var'=tot_`var' if wealth_group==0
	drop tot_`var'
	rename `var' `var'_init
	
	
}

replace prob_EV=. if wealth_group==0
egen m_prob=mean(prob_EV)
replace prob_EV=m_prob if wealth_group==0
drop m_prob
rename prob_EV prob_EV_init

gen subsidy_init=0


**Create the table for the baseline scenario: 

frame change default

rename pv_tax_total tax_pv 

forvalues i=1/4 {
	
frame put _all, into(scen`i')

frame change scen`i'
keep if scenario==`i'

frlink 1:1 wealth_group, frame(init)
frget _all, from(init)

foreach var in co2_ann emb_co2_ann tax_total tax_pv petrol_tax veh_tax {

gen d_`var'=`var'-`var'_init

}

frame put d_co2_ann d_emb_co2_ann d_tax_total d_tax_pv d_petrol_tax d_veh_tax wealth_group CS_change be_change tax_pv petrol_tax veh_tax paid_subsidy, into(Tab10_`i')

frame change Tab10_`i'

replace d_co2_ann=d_co2_ann/1000
replace d_emb_co2_ann=d_emb_co2_ann/1000
gen d_co2_tot=d_co2_ann+d_emb_co2_ann


*Generate experienced CS change: 
gen CSexp_change=CS_change-be_change

local r=0.06
local cap=1/`r' - (1/(`r'*(1+`r')^15))
di `cap'

gen d_wel=(CS_change*1000-(be_change*1000)+d_veh_tax-paid_subsidy+d_tax_pv+(abs(d_co2_tot)*175+d_petrol_tax)*`cap')/1000

gen d_pub=(d_veh_tax-paid_subsidy+d_tax_pv+(d_petrol_tax*`cap'))/1000
gen d_co2_m=abs(d_co2_tot)*`cap'*175/1000

gen TOT_EINK=59514
replace TOT_EINK=87082 if wealth_group==2
replace TOT_EINK=115592 if wealth_group==3
replace TOT_EINK=195182 if wealth_group==4

replace TOT_EINK=114341.3 if wealth_group==0

gen total_tax_pv=-paid_subsidy+veh_tax+tax_pv+(petrol_tax*`cap')

gen total_tax_pv_per_person=total_tax_pv/2307.5
replace total_tax_pv_per_person=total_tax_pv/9230 if wealth_group==0
gen tax_inc=100*total_tax_pv_per_person/(TOT_EINK*`cap')

replace d_petrol_tax=d_petrol_tax/1000
replace d_veh_tax=d_veh_tax/1000
replace d_tax_total=d_tax_total/1000
replace paid_subsidy=paid_subsidy/1000

*Transpose the data to create the desired table: 
gen PanA=.
gen PanB=. 
gen PanC=.
gen PanD=.

keep wealth_group PanA CS_change be_change CSexp_change PanB paid_subsidy d_tax_total d_petrol_tax d_veh_tax tax_inc d_pub PanC d_co2_ann d_emb_co2_ann d_co2_m PanD d_wel

order wealth_group PanA CS_change be_change CSexp_change PanB paid_subsidy d_tax_total d_petrol_tax d_veh_tax tax_inc d_pub PanC d_co2_ann d_emb_co2_ann d_co2_m PanD d_wel 
xpose, clear varname
order _varname, before(v1)

drop in 1

*Create the table: 
mkmat v1 v2 v3 v4 v5 , matrix(Tab10_`i') 

matrix rownames Tab10_`i' = panA cs_d be_d cs_d_tot panB subs tax_d ptax_d veh_tax_d tax_inc d_pub panC co2_d co2_d_emb co2_m panD d_wel
foreach num of numlist 1(1)17 {
foreach numb of numlist 1(1)5 {
matrix Tab10_`i'[`num',`numb']=round(Tab10_`i'[`num',`numb'],.001)
}
}

estadd matrix Tab10_`i', replace
#delim ;
estout matrix(Tab10_`i', fmt(%9.0gc)) using "$root/Results/tables_graphs/welfare_opt_scen`i'.tex", 
style(tex)  mlabels(,none) cells(e(Tab8)(fmt(%9.4f)))
varlabels ( 
panA "\textit{Panel A: Consumer Welfare}"
cs_d "\(\Delta\) Cons. surplus (Decision) (TCHF)"
be_d "\(\Delta\) Belief error (TCHF)"
cs_d_tot "\(\Delta\) Cons. surplus (experienced) (TCHF)"
panB "\textit{Panel B: Public finance}"
subs "Total subsidy (TCHF)"
ptax_d "\(\Delta\) Fuel levy (TCHF p.a)"
tax_d "\(\Delta\) Car registration taxes (TCHF p.a)"
veh_tax_d "\(\Delta\) Vehicle tariffs (TCHF)"
tax_inc "Tax incidence (\%)"
d_pub "\(\Delta\) Public revenue (TCHF)"
panC "\textit{Panel C: Emissions}"
co2_d "\(\Delta\) \(CO_2\) pipe (t p.a.)"
co2_d_emb "\(\Delta\) \(CO_2\) embodied  (t p.a.)"
co2_m "\(\Delta\) \(CO_2\) value (TCHF)"
panD "\textit{Panel D: Overall}"
d_wel "Overall Welfare effect (TCHF)"
             )
collabels ( "Overall"
			"\nth{1} wealth quartile" 
			 "\nth{2} wealth quartile"  
			 "\nth{3} wealth quartile" 
			 "\nth{4} wealth quartile"
			 )
prehead("\begin{tabular}{@{}lccccc@{}} \toprule")  
posthead("\midrule")
prefoot("\bottomrule") 
postfoot("\end{tabular}") replace;
#delim cr


frame change default
frame drop scen`i' 
frame drop Tab10_`i'

}


** Create a graph to compare the different scenarios: 

frame change default
frame put _all, into(graph_stacked)
frame change graph_stacked
keep if wealth_group==0

frlink m:1 wealth_group, frame(init)
frget _all, from(init)

foreach var in co2_ann emb_co2_ann tax_total tax_pv petrol_tax veh_tax {

gen d_`var'=`var'-`var'_init

}

replace d_co2_ann=d_co2_ann/1000
replace d_emb_co2_ann=d_emb_co2_ann/1000
gen d_co2_tot=d_co2_ann+d_emb_co2_ann



gen CSexp_change=CS_change-be_change


local r=0.06
local cap=1/`r' - (1/(`r'*(1+`r')^15))
di `cap'
gen d_pub=(d_veh_tax-paid_subsidy+d_tax_pv+(d_petrol_tax*`cap'))/1000
gen d_co2_m=abs(d_co2_tot)*`cap'*175/1000


keep d_pub d_co2_m CSexp_change scenario 

**Add the data from the subsidy and tax counterfactual: 
frame create subs
frame change subs

import excel "$root/Results/estimation_results/counterfactual/alt2_prob", clear firstrow

gen wealth_group=0

collapse (sum) co2_ann emb_co2_ann tax_total tax_pv petrol_tax veh_tax paid_subsidy, by(wealth_group)

frlink m:1 wealth_group, frame(init)
frget _all, from(init)

foreach var in co2_ann emb_co2_ann tax_total tax_pv petrol_tax veh_tax {

gen d_`var'=`var'-`var'_init

}

replace d_co2_ann=d_co2_ann/1000
replace d_emb_co2_ann=d_emb_co2_ann/1000
gen d_co2_tot=d_co2_ann+d_emb_co2_ann

frame create CS2
frame change CS2

import excel "$root/Results/estimation_results/counterfactual/CS2", clear firstrow

gen CSexp_change=CS_change-be_change

collapse (sum) CSexp_change
gen wealth_group=0

frame change subs 
frlink 1:1 wealth_group, frame(CS2)
frget CSexp_change, from(CS2)

local r=0.06
local cap=1/`r' - (1/(`r'*(1+`r')^15))
di `cap'
gen d_pub=(d_veh_tax-paid_subsidy+d_tax_pv+(d_petrol_tax*`cap'))/1000
gen d_co2_m=abs(d_co2_tot)*`cap'*175/1000

gen scenario=0.2
keep d_pub d_co2_m CSexp_change scenario 


**Add the data from the tax counterfactual: 
frame create tax
frame change tax

import excel "$root/Results/estimation_results/counterfactual/alt3_prob", clear firstrow

gen wealth_group=0
gen paid_subsidy=0
collapse (sum) co2_ann emb_co2_ann tax_total tax_pv petrol_tax veh_tax paid_subsidy, by(wealth_group)

frlink m:1 wealth_group, frame(init)
frget _all, from(init)

foreach var in co2_ann emb_co2_ann tax_total tax_pv petrol_tax veh_tax {

gen d_`var'=`var'-`var'_init

}

replace d_co2_ann=d_co2_ann/1000
replace d_emb_co2_ann=d_emb_co2_ann/1000
gen d_co2_tot=d_co2_ann+d_emb_co2_ann

frame create CS3
frame change CS3

import excel "$root/Results/estimation_results/counterfactual/CS3", clear firstrow

gen CSexp_change=CS_change-be_change

collapse (sum) CSexp_change
gen wealth_group=0

frame change tax 
frlink 1:1 wealth_group, frame(CS3)
frget CSexp_change, from(CS3)

local r=0.06
local cap=1/`r' - (1/(`r'*(1+`r')^15))
di `cap'
gen d_pub=(d_veh_tax-paid_subsidy+d_tax_pv+(d_petrol_tax*`cap'))/1000
gen d_co2_m=abs(d_co2_tot)*`cap'*175/1000

gen scenario=0.1
keep d_pub d_co2_m CSexp_change scenario 


frame change graph_stacked
frameappend subs
frameappend tax


sort scenario

gen d_wel=CSexp_change+d_pub+d_co2_m


replace d_wel=round(d_wel, .01)
replace scenario=scenario+1
replace scenario=0 if _n==1
replace scenario=1 if _n==2
forvalues v=0/5 {
	quietly sum d_wel if scenario==`v'
	local wel_`v'= round(`r(mean)',1)
	
}


set scheme white_tableau
colorpalette Tab10, n(3) nograph gscale

graph hbar d_pub CSexp_change d_co2_m, over(scenario, relabel(1 "Tax Feebate" 2 "EV Subsidy" 3 "Optimal baseline" 4 "Optimal distribution" 5 "Optimal elastic" 6 "Optimal elastic & distribution")) stack legend(pos(6) rows(1) label(1 "Public revenue") label( 3 "CO2 emissions") label( 2 "Consumer Surplus"))  text(1400 95 "Overall: `wel_0' TCHF", placement(e) size(8pt)) text( 1400 76  "Overall: `wel_1' TCHF", placement(e) size(8pt)) text( 1400 58 "Overall: `wel_2' TCHF", placement(e) size(8pt)) text( 1400 42 "Overall: `wel_3' TCHF", placement(e)size(8pt))  text(1400 25 "Overall: `wel_4' TCHF", placement(e) size(8pt)) text(1400 8 "Overall: `wel_5' TCHF", placement(e) size(8pt)) yscale(range(-2000 2200)) ytitle(Monetary value (TCHF)) bar(1, color(`r(p1)')) bar(2, color(`r(p2)')) bar(3, color(`r(p3)'))

*Export the graph: 
save "$root/Results/tables_graphs/main/welfare_stacked_scenarios.png", replace

**Create the additional graphs by wealth group - in total 6 graphs: Share subsidy payments, share tax payments, tax incidence, share in co2 emissions (total), additional EVs, total EVs.

frame change default
frame put _all, into(graph_dist)
frame change graph_dist
drop if wealth_group==0

frame create subs_inc
frame change subs_inc

import excel "$root/Results/estimation_results/counterfactual/alt2_prob_byinc", clear firstrow

gen prob_EV=prob_count2 if drivetype==3

collapse (sum) co2_ann emb_co2_ann petrol_tax paid_subsidy tax_total tax_pv veh_tax prob_EV, by(wealth_group)

gen scenario=.2


frame create tax_inc 
frame change tax_inc


import excel "$root/Results/estimation_results/counterfactual/alt3_prob_byinc", clear firstrow

gen prob_EV=prob_count3 if drivetype==3

gen paid_subsidy=0

collapse (sum) co2_ann emb_co2_ann petrol_tax paid_subsidy tax_total tax_pv veh_tax prob_EV, by(wealth_group)

gen scenario=.1

frame change graph_dist

keep co2_ann emb_co2_ann petrol_tax paid_subsidy tax_total tax_pv veh_tax prob_EV scenario wealth_group 


frameappend subs_inc
frameappend tax_inc


frlink m:1 wealth_group, frame(init)
frget prob_EV_init, from(init)

drop init

gen EV_change=prob_EV-prob_EV_init

gen tot_EV=2307.5*prob_EV
gen add_EV=2307.5*EV_change

gen CO2_tot=co2_ann+ emb_co2_ann
egen tot_emi=sum(CO2_tot), by(scenario)

gen CO2_share=CO2_tot*100/tot_emi

egen tot_subs=sum(paid_subsidy), by(scenario)
gen subs_share=paid_subsidy*100/tot_subs
replace subs_share=0 if missing(paid_subsidy)

local r=0.06
local cap=1/`r' - (1/(`r'*(1+`r')^15))
di `cap'
gen pub_pv=veh_tax+tax_pv+petrol_tax*`cap'
egen tot_pub=sum(pub_pv), by(scenario)
gen share_tax=100*pub_pv/tot_pub

gen TOT_EINK=59514
replace TOT_EINK=87082 if wealth_group==2
replace TOT_EINK=115592 if wealth_group==3
replace TOT_EINK=195182 if wealth_group==4

gen total_tax_pv=-paid_subsidy+veh_tax+tax_pv+(petrol_tax*`cap')

gen total_tax_pv_per_person=total_tax_pv/2307.5
gen tax_inc=100*total_tax_pv_per_person/(TOT_EINK*`cap')

keep wealth_group scenario tax_inc share_tax subs_share CO2_share tot_EV add_EV

replace tot_EV=ceil(tot_EV)
replace add_EV=ceil(add_EV)
reshape wide tax_inc share_tax subs_share CO2_share tot_EV add_EV, i(scenario) j(wealth_group)



	
graph bar tax_inc1 tax_inc2 tax_inc3 tax_inc4, over(scenario, relabel(1 "Tax Feebate" 2 "EV Subsidy" 3 "Optimal baseline" 4 "Optimal distribution" 5 "Optimal elastic" 6 "Optimal elastic & distribution") label(angle(45))) legend(pos(6) rows(2) label(1 "1st wealth quartile") label( 2 "2nd wealth quartile") label( 3 "3rd wealth quartile") label(4 "4th wealth quartile")) ytitle(Tax incidence (%))
	
	
graph export "$root/Results/tables_graphs/appendix/dist_scen_incidence.png", replace
	
	graph bar share_tax1 share_tax2 share_tax3 share_tax4, over(scenario, relabel(1 "Tax Feebate" 2 "EV Subsidy" 3 "Optimal baseline" 4 "Optimal distribution" 5 "Optimal elastic" 6 "Optimal elastic & distribution") label(angle(45))) legend(pos(6) rows(2) label(1 "1st wealth quartile") label( 2 "2nd wealth quartile") label( 3 "3rd wealth quartile") label(4 "4th wealth quartile")) ytitle(Tax share (%))
	
	
graph export "$root/Results/tables_graphs/appendix/dist_scen_taxshare.png", replace

graph bar subs_share1 subs_share2 subs_share3 subs_share4, over(scenario, relabel(1 "Tax Feebate" 2 "EV Subsidy" 3 "Optimal baseline" 4 "Optimal distribution" 5 "Optimal elastic" 6 "Optimal elastic & distribution") label(angle(45))) legend(pos(6) rows(2) label(1 "1st wealth quartile") label( 2 "2nd wealth quartile") label( 3 "3rd wealth quartile") label(4 "4th wealth quartile")) ytitle(Subsidy payments share (%))
	
	
graph export "$root/Results/tables_graphs/appendix/dist_scen_subsshare.png", replace

graph bar CO2_share1 CO2_share2 CO2_share3 CO2_share4 , over(scenario, relabel(1 "Tax Feebate" 2 "EV Subsidy" 3 "Optimal baseline" 4 "Optimal distribution" 5 "Optimal elastic" 6 "Optimal elastic & distribution") label(angle(45))) legend(pos(6) rows(2) label(1 "1st wealth quartile") label( 2 "2nd wealth quartile") label( 3 "3rd wealth quartile") label(4 "4th wealth quartile")) ytitle(CO2 emission share (%))
	
	
graph export "$root/Results/tables_graphs/appendix/dist_scen_co2share.png", replace

graph bar tot_EV1 tot_EV2 tot_EV3 tot_EV4 , over(scenario, relabel(1 "Tax Feebate" 2 "EV Subsidy" 3 "Optimal baseline" 4 "Optimal distribution" 5 "Optimal elastic" 6 "Optimal elastic & distribution") label(angle(45))) legend(pos(6) rows(2) label(1 "1st wealth quartile") label( 2 "2nd wealth quartile") label( 3 "3rd wealth quartile") label(4 "4th wealth quartile")) ytitle(Total EVs)
	
	
graph export "$root/Results/tables_graphs/appendix/dist_scen_evtot.png", replace

graph bar add_EV1 add_EV2 add_EV3 add_EV4 , over(scenario, relabel(1 "Tax Feebate" 2 "EV Subsidy" 3 "Optimal baseline" 4 "Optimal distribution" 5 "Optimal elastic" 6 "Optimal elastic & distribution") label(angle(45))) legend(pos(6) rows(2) label(1 "1st wealth quartile") label( 2 "2nd wealth quartile") label( 3 "3rd wealth quartile") label(4 "4th wealth quartile")) ytitle(Additional EVs)
	
	
graph export "$root/Results/tables_graphs/appendix/dist_scen_evadd.png", replace
